clc;
clear;
close all;
rng(222);

% 遗传算法参数
popSize = 100;     % 种群规模
nGenerations = 20; % 代数
crossoverRate = 0.9; % 交叉概率
mutationRate = 0.1;  % 变异概率
n = popSize;        % 个体数量
m = 2;              % 目标数

% 初始化种群
pop = randn(n,m);

for generation = 1:nGenerations
    % 计算目标函数值
    obj = mubi(pop(:,1),pop(:,2)); 
    
    % 计算非支配排序
    [F, ndx] = fast_nondominated_sort(obj);
    
    % 计算拥挤距离
    crowdingDist = calculateCrowdingDistance(obj, F);
    
    % 选择
    selectedPop = selection(pop, F, crowdingDist, popSize);
    
    % 交叉
    offspring = crossover(selectedPop, crossoverRate);
    
    % 变异
    offspring = mutation(offspring, mutationRate);
    
    % 更新种群
    pop = updatePopulation(pop, offspring, popSize);
    
    % 打印当前种群的分布
    figure(2);
    scatter(pop(:,1), pop(:,2), 'filled');
    title(['Population at Generation ', num2str(generation)]);
    xlabel('x');
    ylabel('y');
    pause(0.1); % 暂停以便查看图形
end

% 绘制最终的 Pareto 前沿
figure;
hold on;
obj = mubi(pop(:,1), pop(:,2)); % 计算目标函数值
[F, ~] = fast_nondominated_sort(obj); % 重新计算前沿
colors = lines(length(F)); % 为不同层次分配不同颜色
for i = 1:length(F)
    front_points = pop(F{i}, :);
    [~, sort_idx] = sort(front_points(:, 1).^2); % 根据目标排序
    sorted_front_points = front_points(sort_idx, :);
    plot(sorted_front_points(:, 1), sorted_front_points(:, 2), 'o-', 'Color', colors(i, :), 'DisplayName', sprintf('Front %d', i));
end
xlabel('x');
ylabel('y');
title('Pareto Front');
grid on;
legend show;

%% 目标函数定义
function obj=mubi(x,y)
    obj(:,1) = y.^2+x.^2;     % min x^2
    obj(:,2) = -y.^2;    % max y^2 (注意：这是通过取负来转换为最小化问题)
end

%% 非支配排序函数
function [F, ndx] = fast_nondominated_sort(obj)
    [n, m] = size(obj);
    F = cell(1, n);
    S = cell(1, n);
    nPoint = zeros(1, n);
    ndx = zeros(1, n); % 初始化 ndx

    for i = 1:n
        S{i} = [];
        nPoint(i) = 0;
        for j = 1:n
            if istro(obj(i, :), obj(j, :))
                S{i} = [S{i}, j];
            elseif istro(obj(j, :), obj(i, :))
                nPoint(i) = nPoint(i) + 1;
            end
        end
        if nPoint(i) == 0
            F{1} = [F{1}, i];
        end
    end

    ifs = 1;
    while ~isempty(F{ifs})
        Q = [];
        for i = F{ifs}
            for j = S{i}
                nPoint(j) = nPoint(j) - 1;
                if nPoint(j) == 0
                    Q = [Q, j];
                    ndx(j) = ifs + 1;
                end
            end
        end
        ifs = ifs + 1;
        F{ifs} = Q;
    end
    F = F(~cellfun('isempty', F)); % 去除空的 cell
end

%% 拥挤距离计算
function crowdDist = calculateCrowdingDistance(obj, F)
    [n, m] = size(obj);
    crowdDist = zeros(n, 1);
    
    for i = 1:length(F)
        front = F{i};  % 当前前沿的个体索引
        frontSize = length(front);
        
        if frontSize == 0
            continue;
        end
        
        dist = zeros(frontSize, 1);  % 初始化当前前沿的拥挤距离
        for j = 1:m
            % 按第 j 个目标值排序
            [~, sortIdx] = sort(obj(front, j));
            sortedFront = front(sortIdx);
            
            % 设置边界个体的拥挤距离为无穷大
            dist(1) = Inf;
            dist(end) = Inf;
            
            % 计算中间个体的拥挤距离
            for k = 2:frontSize-1
                dist(k) = dist(k) + (obj(sortedFront(k+1), j) - obj(sortedFront(k-1), j)) / (max(obj(:, j)) - min(obj(:, j)));
            end
        end
        
        % 将拥挤距离赋给当前前沿的个体
        crowdDist(front) = dist;
    end
end

%% 判断支配关系
function flag = istro(x, y)
    flag = all(x <= y) && any(x < y);
end

%% 选择操作（基于非支配排序和拥挤距离）
function selectedPop = selection(pop, F, crowdingDist, popSize)
    selectedPop = [];
    i = 1;
    while length(selectedPop) + length(F{i}) < popSize
        selectedPop = [selectedPop; pop(F{i}, :)];  % 选择整个前沿的个体
        i = i + 1;
    end
    
    % 最后一层按拥挤距离排序，选择剩余个体
    remaining = popSize - length(selectedPop);
    if remaining > 0
        lastFront = F{i};
        [~, sortIdx] = sort(crowdingDist(lastFront), 'descend');  % 按拥挤距离排序
        selected = lastFront(sortIdx(1:remaining));  % 选择剩余的个体
        selectedPop = [selectedPop; pop(selected, :)];
    end
end

%% 交叉操作
function offspring = crossover(pop, crossoverRate)
    [n, m] = size(pop);
    numCrossovers = round(crossoverRate * n / 2);
    offspring = pop;
    
    for i = 1:numCrossovers
        parents = randperm(n, 2);
        p1 = pop(parents(1), :);
        p2 = pop(parents(2), :);
        point = randi(m-1); % 交叉点
        offspring(i*2-1, :) = [p1(1:point), p2(point+1:end)];
        offspring(i*2, :) = [p2(1:point), p1(point+1:end)];
    end
end

%% 变异操作
function mutatedPop = mutation(pop, mutationRate)
    [n, m] = size(pop);
    mutatedPop = pop;
    numMutations = round(mutationRate * n * m);
    
    for i = 1:numMutations
        individual = randi(n);
        gene = randi(m);
        mutatedPop(individual, gene) = rand() * 20 - 10; % 保持在 [-10, 10] 范围内
    end
end

%% 更新种群
function newPop = updatePopulation(oldPop, offspring, popSize)
    combinedPop = [oldPop; offspring];
    obj = mubi(combinedPop(:,1), combinedPop(:,2)); % 计算目标函数值
    [F, ndx] = fast_nondominated_sort(obj);
    crowdingDist = calculateCrowdingDistance(obj, F);
    newPop = [];
    i = 1;
    while length(newPop) + length(F{i}) <= popSize
        newPop = [newPop; combinedPop(F{i}, :)];  % 添加整个前沿
        i = i + 1;
    end
    
    % 对最后一个前沿，按拥挤距离排序
    remaining = popSize - length(newPop);  % 还需要补充的个体数量
    if remaining > 0
        lastFront = F{i};  % 最后一层
        [~, sortIdx] = sort(crowdingDist(lastFront), 'descend');  % 拥挤距离从大到小排序
        selected = lastFront(sortIdx(1:remaining));  % 选择前 remaining 个个体
        newPop = [newPop; combinedPop(selected, :)];  % 加入到新种群中
    end
end
